home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / coulomb.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  39.5 KB  |  1,416 lines
  1. /* specfunc/coulomb.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. /* Evaluation of Coulomb wave functions F_L(eta, x), G_L(eta, x),
  23.  * and their derivatives. A combination of Steed's method, asymptotic
  24.  * results, and power series.
  25.  *
  26.  * Steed's method:
  27.  *  [Barnett, CPC 21, 297 (1981)]
  28.  * Power series and other methods:
  29.  *  [Biedenharn et al., PR 97, 542 (1954)]
  30.  *  [Bardin et al., CPC 3, 73 (1972)]
  31.  *  [Abad+Sesma, CPC 71, 110 (1992)]
  32.  */
  33. #include <config.h>
  34. #include <gsl/gsl_math.h>
  35. #include <gsl/gsl_errno.h>
  36. #include <gsl/gsl_sf_exp.h>
  37. #include <gsl/gsl_sf_psi.h>
  38. #include <gsl/gsl_sf_airy.h>
  39. #include <gsl/gsl_sf_pow_int.h>
  40. #include <gsl/gsl_sf_gamma.h>
  41. #include <gsl/gsl_sf_coulomb.h>
  42.  
  43. #include "error.h"
  44.  
  45. /* the L=0 normalization constant
  46.  * [Abramowitz+Stegun 14.1.8]
  47.  */
  48. static
  49. double
  50. C0sq(double eta)
  51. {
  52.   double twopieta = 2.0*M_PI*eta;
  53.  
  54.   if(fabs(eta) < GSL_DBL_EPSILON) {
  55.     return 1.0;
  56.   }
  57.   else if(twopieta > GSL_LOG_DBL_MAX) {
  58.     return 0.0;
  59.   }
  60.   else {
  61.     gsl_sf_result scale;
  62.     gsl_sf_expm1_e(twopieta, &scale);
  63.     return twopieta/scale.val;
  64.   }
  65. }
  66.  
  67.  
  68. /* the full definition of C_L(eta) for any valid L and eta
  69.  * [Abramowitz and Stegun 14.1.7]
  70.  * This depends on the complex gamma function. For large
  71.  * arguments the phase of the complex gamma function is not
  72.  * very accurately determined. However the modulus is, and that
  73.  * is all that we need to calculate C_L.
  74.  *
  75.  * This is not valid for L <= -3/2  or  L = -1.
  76.  */
  77. static
  78. int
  79. CLeta(double L, double eta, gsl_sf_result * result)
  80. {
  81.   gsl_sf_result ln1; /* log of numerator Gamma function */
  82.   gsl_sf_result ln2; /* log of denominator Gamma function */
  83.   double sgn = 1.0;
  84.   double arg_val, arg_err;
  85.  
  86.   if(fabs(eta/(L+1.0)) < GSL_DBL_EPSILON) {
  87.     gsl_sf_lngamma_e(L+1.0, &ln1);
  88.   }
  89.   else {
  90.     gsl_sf_result p1;                 /* phase of numerator Gamma -- not used */
  91.     gsl_sf_lngamma_complex_e(L+1.0, eta, &ln1, &p1); /* should be ok */
  92.   }
  93.  
  94.   gsl_sf_lngamma_e(2.0*(L+1.0), &ln2);
  95.   if(L < -1.0) sgn = -sgn;
  96.  
  97.   arg_val  = L*M_LN2 - 0.5*eta*M_PI + ln1.val - ln2.val;
  98.   arg_err  = ln1.err + ln2.err;
  99.   arg_err += GSL_DBL_EPSILON * (fabs(L*M_LN2) + fabs(0.5*eta*M_PI));
  100.   return gsl_sf_exp_err_e(arg_val, arg_err, result);
  101. }
  102.  
  103.  
  104. int
  105. gsl_sf_coulomb_CL_e(double lam, double eta, gsl_sf_result * result)
  106. {
  107.   /* CHECK_POINTER(result) */
  108.  
  109.   if(lam <= -1.0) {
  110.     DOMAIN_ERROR(result);
  111.   }
  112.   else if(fabs(lam) < GSL_DBL_EPSILON) {
  113.     /* saves a calculation of complex_lngamma(), otherwise not necessary */
  114.     result->val = sqrt(C0sq(eta));
  115.     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  116.     return GSL_SUCCESS;
  117.   }
  118.   else {
  119.     return CLeta(lam, eta, result);
  120.   }
  121. }
  122.  
  123.  
  124. /* cl[0] .. cl[kmax] = C_{lam_min}(eta) .. C_{lam_min+kmax}(eta)
  125.  */
  126. int
  127. gsl_sf_coulomb_CL_array(double lam_min, int kmax, double eta, double * cl)
  128. {
  129.   int k;
  130.   gsl_sf_result cl_0;
  131.   gsl_sf_coulomb_CL_e(lam_min, eta, &cl_0);
  132.   cl[0] = cl_0.val;
  133.  
  134.   for(k=1; k<=kmax; k++) {
  135.     double L = lam_min + k;
  136.     cl[k] = cl[k-1] * sqrt(L*L + eta*eta)/(L*(2.0*L+1.0));
  137.   }
  138.  
  139.   return GSL_SUCCESS;
  140. }
  141.  
  142.  
  143. /* Evaluate the series for Phi_L(eta,x) and Phi_L*(eta,x)
  144.  * [Abramowitz+Stegun 14.1.5]
  145.  * [Abramowitz+Stegun 14.1.13]
  146.  *
  147.  * The sequence of coefficients A_k^L is
  148.  * manifestly well-controlled for L >= -1/2
  149.  * and eta < 10.
  150.  *
  151.  * This makes sense since this is the region
  152.  * away from threshold, and you expect
  153.  * the evaluation to become easier as you
  154.  * get farther from threshold.
  155.  *
  156.  * Empirically, this is quite well-behaved for
  157.  *   L >= -1/2
  158.  *   eta < 10
  159.  *   x   < 10
  160.  */
  161. #if 0
  162. static
  163. int
  164. coulomb_Phi_series(const double lam, const double eta, const double x,
  165.                    double * result, double * result_star)
  166. {
  167.   int kmin =   5;
  168.   int kmax = 200;
  169.   int k;
  170.   double Akm2 = 1.0;
  171.   double Akm1 = eta/(lam+1.0);
  172.   double Ak;
  173.  
  174.   double xpow = x;
  175.   double sum  = Akm2 + Akm1*x;
  176.   double sump = (lam+1.0)*Akm2 + (lam+2.0)*Akm1*x;
  177.   double prev_abs_del   = fabs(Akm1*x);
  178.   double prev_abs_del_p = (lam+2.0) * prev_abs_del;
  179.  
  180.   for(k=2; k<kmax; k++) {
  181.     double del;
  182.     double del_p;
  183.     double abs_del;
  184.     double abs_del_p;
  185.  
  186.     Ak = (2.0*eta*Akm1 - Akm2)/(k*(2.0*lam + 1.0 + k));
  187.  
  188.     xpow *= x;
  189.     del   = Ak*xpow;
  190.     del_p = (k+lam+1.0)*del;
  191.     sum  += del;
  192.     sump += del_p;
  193.  
  194.     abs_del   = fabs(del);
  195.     abs_del_p = fabs(del_p);
  196.  
  197.     if(          abs_del/(fabs(sum)+abs_del)          < GSL_DBL_EPSILON
  198.        &&   prev_abs_del/(fabs(sum)+prev_abs_del)     < GSL_DBL_EPSILON
  199.        &&      abs_del_p/(fabs(sump)+abs_del_p)       < GSL_DBL_EPSILON
  200.        && prev_abs_del_p/(fabs(sump)+prev_abs_del_p)  < GSL_DBL_EPSILON
  201.        && k > kmin
  202.        ) break;
  203.  
  204.     /* We need to keep track of the previous delta because when
  205.      * eta is near zero the odd terms of the sum are very small
  206.      * and this could lead to premature termination.
  207.      */
  208.     prev_abs_del   = abs_del;
  209.     prev_abs_del_p = abs_del_p;
  210.  
  211.     Akm2 = Akm1;
  212.     Akm1 = Ak;
  213.   }
  214.  
  215.   *result      = sum;
  216.   *result_star = sump;
  217.  
  218.   if(k==kmax) {
  219.     GSL_ERROR ("error", GSL_EMAXITER);
  220.   }
  221.   else {
  222.     return GSL_SUCCESS;
  223.   }
  224. }
  225. #endif /* 0 */
  226.  
  227.  
  228. /* Determine the connection phase, phi_lambda.
  229.  * See coulomb_FG_series() below. We have
  230.  * to be careful about sin(phi)->0. Note that
  231.  * there is an underflow condition for large 
  232.  * positive eta in any case.
  233.  */
  234. static
  235. int
  236. coulomb_connection(const double lam, const double eta,
  237.                    double * cos_phi, double * sin_phi)
  238. {
  239.   if(eta > -GSL_LOG_DBL_MIN/2.0*M_PI-1.0) {
  240.     *cos_phi = 1.0;
  241.     *sin_phi = 0.0;
  242.     GSL_ERROR ("error", GSL_EUNDRFLW);
  243.   }
  244.   else if(eta > -GSL_LOG_DBL_EPSILON/(4.0*M_PI)) {
  245.     const double eps = 2.0 * exp(-2.0*M_PI*eta);
  246.     const double tpl = tan(M_PI * lam);
  247.     const double dth = eps * tpl / (tpl*tpl + 1.0);
  248.     *cos_phi = -1.0 + 0.5 * dth*dth;
  249.     *sin_phi = -dth;
  250.     return GSL_SUCCESS;
  251.   }
  252.   else {
  253.     double X   = tanh(M_PI * eta) / tan(M_PI * lam);
  254.     double phi = -atan(X) - (lam + 0.5) * M_PI;
  255.     *cos_phi = cos(phi);
  256.     *sin_phi = sin(phi);
  257.     return GSL_SUCCESS;
  258.   }
  259. }
  260.  
  261.  
  262. /* Evaluate the Frobenius series for F_lam(eta,x) and G_lam(eta,x).
  263.  * Homegrown algebra. Evaluates the series for F_{lam} and
  264.  * F_{-lam-1}, then uses
  265.  *    G_{lam} = (F_{lam} cos(phi) - F_{-lam-1}) / sin(phi)
  266.  * where
  267.  *    phi = Arg[Gamma[1+lam+I eta]] - Arg[Gamma[-lam + I eta]] - (lam+1/2)Pi
  268.  *        = Arg[Sin[Pi(-lam+I eta)] - (lam+1/2)Pi
  269.  *        = atan2(-cos(lam Pi)sinh(eta Pi), -sin(lam Pi)cosh(eta Pi)) - (lam+1/2)Pi
  270.  *
  271.  *        = -atan(X) - (lam+1/2) Pi,  X = tanh(eta Pi)/tan(lam Pi)
  272.  *
  273.  * Not appropriate for lam <= -1/2, lam = 0, or lam >= 1/2.
  274.  */
  275. static
  276. int
  277. coulomb_FG_series(const double lam, const double eta, const double x,
  278.                   gsl_sf_result * F, gsl_sf_result * G)
  279. {
  280.   const int max_iter = 800;
  281.   gsl_sf_result ClamA;
  282.   gsl_sf_result ClamB;
  283.   int stat_A = CLeta(lam, eta, &ClamA);
  284.   int stat_B = CLeta(-lam-1.0, eta, &ClamB);
  285.   const double tlp1 = 2.0*lam + 1.0;
  286.   const double pow_x = pow(x, lam);
  287.   double cos_phi_lam;
  288.   double sin_phi_lam;
  289.  
  290.   double uA_mm2 = 1.0;                  /* uA sum is for F_{lam} */
  291.   double uA_mm1 = x*eta/(lam+1.0);
  292.   double uA_m;
  293.   double uB_mm2 = 1.0;                  /* uB sum is for F_{-lam-1} */
  294.   double uB_mm1 = -x*eta/lam;
  295.   double uB_m;
  296.   double A_sum = uA_mm2 + uA_mm1;
  297.   double B_sum = uB_mm2 + uB_mm1;
  298.   double A_abs_del_prev = fabs(A_sum);
  299.   double B_abs_del_prev = fabs(B_sum);
  300.   gsl_sf_result FA, FB;
  301.   int m = 2;
  302.  
  303.   int stat_conn = coulomb_connection(lam, eta, &cos_phi_lam, &sin_phi_lam);
  304.  
  305.   if(stat_conn == GSL_EUNDRFLW) {
  306.     F->val = 0.0;  /* FIXME: should this be set to Inf too like G? */
  307.     F->err = 0.0;
  308.     OVERFLOW_ERROR(G);
  309.   }
  310.  
  311.   while(m < max_iter) {
  312.     double abs_dA;
  313.     double abs_dB;
  314.     uA_m = x*(2.0*eta*uA_mm1 - x*uA_mm2)/(m*(m+tlp1));
  315.     uB_m = x*(2.0*eta*uB_mm1 - x*uB_mm2)/(m*(m-tlp1));
  316.     A_sum += uA_m;
  317.     B_sum += uB_m;
  318.     abs_dA = fabs(uA_m);
  319.     abs_dB = fabs(uB_m);
  320.     if(m > 15) {
  321.       /* Don't bother checking until we have gone out a little ways;
  322.        * a minor optimization. Also make sure to check both the
  323.        * current and the previous increment because the odd and even
  324.        * terms of the sum can have very different behaviour, depending
  325.        * on the value of eta.
  326.        */
  327.       double max_abs_dA = GSL_MAX(abs_dA, A_abs_del_prev);
  328.       double max_abs_dB = GSL_MAX(abs_dB, B_abs_del_prev);
  329.       double abs_A = fabs(A_sum);
  330.       double abs_B = fabs(B_sum);
  331.       if(   max_abs_dA/(max_abs_dA + abs_A) < 4.0*GSL_DBL_EPSILON
  332.          && max_abs_dB/(max_abs_dB + abs_B) < 4.0*GSL_DBL_EPSILON
  333.          ) break;
  334.     }
  335.     A_abs_del_prev = abs_dA;
  336.     B_abs_del_prev = abs_dB;
  337.     uA_mm2 = uA_mm1;
  338.     uA_mm1 = uA_m;
  339.     uB_mm2 = uB_mm1;
  340.     uB_mm1 = uB_m;
  341.     m++;
  342.   }
  343.  
  344.   FA.val = A_sum * ClamA.val * pow_x * x;
  345.   FA.err = fabs(A_sum) * ClamA.err * pow_x * x + 2.0*GSL_DBL_EPSILON*fabs(FA.val);
  346.   FB.val = B_sum * ClamB.val / pow_x;
  347.   FB.err = fabs(B_sum) * ClamB.err / pow_x + 2.0*GSL_DBL_EPSILON*fabs(FB.val);
  348.  
  349.   F->val = FA.val;
  350.   F->err = FA.err;
  351.  
  352.   G->val = (FA.val * cos_phi_lam - FB.val)/sin_phi_lam;
  353.   G->err = (FA.err * fabs(cos_phi_lam) + FB.err)/fabs(sin_phi_lam);
  354.  
  355.   if(m >= max_iter)
  356.     GSL_ERROR ("error", GSL_EMAXITER);
  357.   else
  358.     return GSL_ERROR_SELECT_2(stat_A, stat_B);
  359. }
  360.  
  361.  
  362. /* Evaluate the Frobenius series for F_0(eta,x) and G_0(eta,x).
  363.  * See [Bardin et al., CPC 3, 73 (1972), (14)-(17)];
  364.  * note the misprint in (17): nu_0=1 is correct, not nu_0=0.
  365.  */
  366. static
  367. int
  368. coulomb_FG0_series(const double eta, const double x,
  369.                    gsl_sf_result * F, gsl_sf_result * G)
  370. {
  371.   const int max_iter = 800;
  372.   const double x2  = x*x;
  373.   const double tex = 2.0*eta*x;
  374.   gsl_sf_result C0;
  375.   int stat_CL = CLeta(0.0, eta, &C0);
  376.   gsl_sf_result r1pie;
  377.   int psi_stat = gsl_sf_psi_1piy_e(eta, &r1pie);
  378.   double u_mm2 = 0.0;  /* u_0 */
  379.   double u_mm1 = x;    /* u_1 */
  380.   double u_m;
  381.   double v_mm2 = 1.0;                               /* nu_0 */
  382.   double v_mm1 = tex*(2.0*M_EULER-1.0+r1pie.val);   /* nu_1 */
  383.   double v_m;
  384.   double u_sum = u_mm2 + u_mm1;
  385.   double v_sum = v_mm2 + v_mm1;
  386.   double u_abs_del_prev = fabs(u_sum);
  387.   double v_abs_del_prev = fabs(v_sum);
  388.   int m = 2;
  389.   double u_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(u_sum);
  390.   double v_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(v_sum);
  391.   double ln2x = log(2.0*x);
  392.  
  393.   while(m < max_iter) {
  394.     double abs_du;
  395.     double abs_dv;
  396.     double m_mm1 = m*(m-1.0);
  397.     u_m = (tex*u_mm1 - x2*u_mm2)/m_mm1;
  398.     v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*eta*(2*m-1)*u_m)/m_mm1;
  399.     u_sum += u_m;
  400.     v_sum += v_m;
  401.     abs_du = fabs(u_m);
  402.     abs_dv = fabs(v_m);
  403.     u_sum_err += 2.0 * GSL_DBL_EPSILON * abs_du;
  404.     v_sum_err += 2.0 * GSL_DBL_EPSILON * abs_dv;
  405.     if(m > 15) {
  406.       /* Don't bother checking until we have gone out a little ways;
  407.        * a minor optimization. Also make sure to check both the
  408.        * current and the previous increment because the odd and even
  409.        * terms of the sum can have very different behaviour, depending
  410.        * on the value of eta.
  411.        */
  412.       double max_abs_du = GSL_MAX(abs_du, u_abs_del_prev);
  413.       double max_abs_dv = GSL_MAX(abs_dv, v_abs_del_prev);
  414.       double abs_u = fabs(u_sum);
  415.       double abs_v = fabs(v_sum);
  416.       if(   max_abs_du/(max_abs_du + abs_u) < 40.0*GSL_DBL_EPSILON
  417.          && max_abs_dv/(max_abs_dv + abs_v) < 40.0*GSL_DBL_EPSILON
  418.          ) break;
  419.     }
  420.     u_abs_del_prev = abs_du;
  421.     v_abs_del_prev = abs_dv;
  422.     u_mm2 = u_mm1;
  423.     u_mm1 = u_m;
  424.     v_mm2 = v_mm1;
  425.     v_mm1 = v_m;
  426.     m++;
  427.   }
  428.  
  429.   F->val  = C0.val * u_sum;
  430.   F->err  = C0.err * fabs(u_sum);
  431.   F->err += fabs(C0.val) * u_sum_err;
  432.   F->err += 2.0 * GSL_DBL_EPSILON * fabs(F->val);
  433.  
  434.   G->val  = (v_sum + 2.0*eta*u_sum * ln2x) / C0.val;
  435.   G->err  = (fabs(v_sum) + fabs(2.0*eta*u_sum * ln2x)) / fabs(C0.val) * fabs(C0.err/C0.val);
  436.   G->err += (v_sum_err + fabs(2.0*eta*u_sum_err*ln2x)) / fabs(C0.val);
  437.   G->err += 2.0 * GSL_DBL_EPSILON * fabs(G->val);
  438.  
  439.   if(m == max_iter)
  440.     GSL_ERROR ("error", GSL_EMAXITER);
  441.   else
  442.     return GSL_ERROR_SELECT_2(psi_stat, stat_CL);
  443. }
  444.  
  445.  
  446. /* Evaluate the Frobenius series for F_{-1/2}(eta,x) and G_{-1/2}(eta,x).
  447.  * Homegrown algebra.
  448.  */
  449. static
  450. int
  451. coulomb_FGmhalf_series(const double eta, const double x,
  452.                        gsl_sf_result * F, gsl_sf_result * G)
  453. {
  454.   const int max_iter = 800;
  455.   const double rx  = sqrt(x);
  456.   const double x2  = x*x;
  457.   const double tex = 2.0*eta*x;
  458.   gsl_sf_result Cmhalf;
  459.   int stat_CL = CLeta(-0.5, eta, &Cmhalf);
  460.   double u_mm2 = 1.0;                      /* u_0 */
  461.   double u_mm1 = tex * u_mm2;              /* u_1 */
  462.   double u_m;
  463.   double v_mm2, v_mm1, v_m;
  464.   double f_sum, g_sum;
  465.   double tmp1;
  466.   gsl_sf_result rpsi_1pe;
  467.   gsl_sf_result rpsi_1p2e;
  468.   int m = 2;
  469.  
  470.   gsl_sf_psi_1piy_e(eta,     &rpsi_1pe);
  471.   gsl_sf_psi_1piy_e(2.0*eta, &rpsi_1p2e);
  472.  
  473.   v_mm2 = 2.0*M_EULER - M_LN2 - rpsi_1pe.val + 2.0*rpsi_1p2e.val;
  474.   v_mm1 = tex*(v_mm2 - 2.0*u_mm2);
  475.  
  476.   f_sum = u_mm2 + u_mm1;
  477.   g_sum = v_mm2 + v_mm1;
  478.  
  479.   while(m < max_iter) {
  480.     double m2 = m*m;
  481.     u_m = (tex*u_mm1 - x2*u_mm2)/m2;
  482.     v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*m*u_m)/m2;
  483.     f_sum += u_m;
  484.     g_sum += v_m;
  485.     if(   f_sum != 0.0
  486.        && g_sum != 0.0
  487.        && (fabs(u_m/f_sum) + fabs(v_m/g_sum) < 10.0*GSL_DBL_EPSILON)) break;
  488.     u_mm2 = u_mm1;
  489.     u_mm1 = u_m;
  490.     v_mm2 = v_mm1;
  491.     v_mm1 = v_m;
  492.     m++;
  493.   }
  494.   
  495.   F->val = Cmhalf.val * rx * f_sum;
  496.   F->err = Cmhalf.err * fabs(rx * f_sum) + 2.0*GSL_DBL_EPSILON*fabs(F->val);
  497.  
  498.   tmp1 = f_sum*log(x);
  499.   G->val = -rx*(tmp1 + g_sum)/Cmhalf.val;
  500.   G->err = fabs(rx)*(fabs(tmp1) + fabs(g_sum))/fabs(Cmhalf.val) * fabs(Cmhalf.err/Cmhalf.val);
  501.  
  502.   if(m == max_iter)
  503.     GSL_ERROR ("error", GSL_EMAXITER);
  504.   else
  505.     return stat_CL;
  506. }
  507.  
  508.  
  509. /* Evolve the backwards recurrence for F,F'.
  510.  *
  511.  *    F_{lam-1}  = (S_lam F_lam + F_lam') / R_lam
  512.  *    F_{lam-1}' = (S_lam F_{lam-1} - R_lam F_lam)
  513.  * where
  514.  *    R_lam = sqrt(1 + (eta/lam)^2)
  515.  *    S_lam = lam/x + eta/lam
  516.  *
  517.  */
  518. static
  519. int
  520. coulomb_F_recur(double lam_min, int kmax,
  521.                 double eta, double x,
  522.                 double F_lam_max, double Fp_lam_max,
  523.         double * F_lam_min, double * Fp_lam_min
  524.                 )
  525. {
  526.   double x_inv = 1.0/x;
  527.   double fcl = F_lam_max;
  528.   double fpl = Fp_lam_max;
  529.   double lam_max = lam_min + kmax;
  530.   double lam = lam_max;
  531.   int k;
  532.  
  533.   for(k=kmax-1; k>=0; k--) {
  534.     double el = eta/lam;
  535.     double rl = sqrt(1.0 + el*el);
  536.     double sl = el  + lam*x_inv;
  537.     double fc_lm1;
  538.     fc_lm1 = (fcl*sl + fpl)/rl;
  539.     fpl    =  fc_lm1*sl - fcl*rl;
  540.     fcl    =  fc_lm1;
  541.     lam -= 1.0;
  542.   }
  543.  
  544.   *F_lam_min  = fcl;
  545.   *Fp_lam_min = fpl;  
  546.   return GSL_SUCCESS;
  547. }
  548.  
  549.  
  550. /* Evolve the forward recurrence for G,G'.
  551.  *
  552.  *   G_{lam+1}  = (S_lam G_lam - G_lam')/R_lam
  553.  *   G_{lam+1}' = R_{lam+1} G_lam - S_lam G_{lam+1}
  554.  *
  555.  * where S_lam and R_lam are as above in the F recursion.
  556.  */
  557. static
  558. int
  559. coulomb_G_recur(const double lam_min, const int kmax,
  560.                 const double eta, const double x,
  561.                 const double G_lam_min, const double Gp_lam_min,
  562.         double * G_lam_max, double * Gp_lam_max
  563.                 )
  564. {
  565.   double x_inv = 1.0/x;
  566.   double gcl = G_lam_min;
  567.   double gpl = Gp_lam_min;
  568.   double lam = lam_min + 1.0;
  569.   int k;
  570.  
  571.   for(k=1; k<=kmax; k++) {
  572.     double el = eta/lam;
  573.     double rl = sqrt(1.0 + el*el);
  574.     double sl = el + lam*x_inv;
  575.     double gcl1 = (sl*gcl - gpl)/rl;
  576.     gpl   = rl*gcl - sl*gcl1;
  577.     gcl   = gcl1;
  578.     lam += 1.0;
  579.   }
  580.   
  581.   *G_lam_max  = gcl;
  582.   *Gp_lam_max = gpl;
  583.   return GSL_SUCCESS;
  584. }
  585.  
  586.  
  587. /* Evaluate the first continued fraction, giving
  588.  * the ratio F'/F at the upper lambda value.
  589.  * We also determine the sign of F at that point,
  590.  * since it is the sign of the last denominator
  591.  * in the continued fraction.
  592.  */
  593. static
  594. int
  595. coulomb_CF1(double lambda,
  596.             double eta, double x,
  597.             double * fcl_sign,
  598.         double * result,
  599.         int * count
  600.             )
  601. {
  602.   const double CF1_small = 1.e-30;
  603.   const double CF1_abort = 1.0e+05;
  604.   const double CF1_acc   = 2.0*GSL_DBL_EPSILON;
  605.   const double x_inv     = 1.0/x;
  606.   const double px        = lambda + 1.0 + CF1_abort;
  607.  
  608.   double pk = lambda + 1.0;
  609.   double F  = eta/pk + pk*x_inv;
  610.   double D, C;
  611.   double df;
  612.  
  613.   *fcl_sign = 1.0;
  614.   *count = 0;
  615.  
  616.   if(fabs(F) < CF1_small) F = CF1_small;
  617.   D = 0.0;
  618.   C = F;
  619.  
  620.   do {
  621.     double pk1 = pk + 1.0;
  622.     double ek  = eta / pk;
  623.     double rk2 = 1.0 + ek*ek;
  624.     double tk  = (pk + pk1)*(x_inv + ek/pk1);
  625.     D   =  tk - rk2 * D;
  626.     C   =  tk - rk2 / C;
  627.     if(fabs(C) < CF1_small) C = CF1_small;
  628.     if(fabs(D) < CF1_small) D = CF1_small;
  629.     D = 1.0/D;
  630.     df = D * C;
  631.     F  = F * df;
  632.     if(D < 0.0) {
  633.       /* sign of result depends on sign of denominator */
  634.       *fcl_sign = - *fcl_sign;
  635.     }
  636.     pk = pk1;
  637.     if( pk > px ) {
  638.       *result = F;
  639.       GSL_ERROR ("error", GSL_ERUNAWAY);
  640.     }
  641.     ++(*count);
  642.   }
  643.   while(fabs(df-1.0) > CF1_acc);
  644.   
  645.   *result = F;
  646.   return GSL_SUCCESS;
  647. }
  648.  
  649.  
  650. #if 0
  651. static
  652. int
  653. old_coulomb_CF1(const double lambda,
  654.                 double eta, double x,
  655.                 double * fcl_sign,
  656.             double * result
  657.                 )
  658. {
  659.   const double CF1_abort = 1.e5;
  660.   const double CF1_acc   = 10.0*GSL_DBL_EPSILON;
  661.   const double x_inv     = 1.0/x;
  662.   const double px        = lambda + 1.0 + CF1_abort;
  663.   
  664.   double pk = lambda + 1.0;
  665.   
  666.   double D;
  667.   double df;
  668.  
  669.   double F;
  670.   double p;
  671.   double pk1;
  672.   double ek;
  673.   
  674.   double fcl = 1.0;
  675.  
  676.   double tk;
  677.  
  678.   while(1) {
  679.     ek = eta/pk;
  680.     F = (ek + pk*x_inv)*fcl + (fcl - 1.0)*x_inv;
  681.     pk1 = pk + 1.0;
  682.     if(fabs(eta*x + pk*pk1) > CF1_acc) break;
  683.     fcl = (1.0 + ek*ek)/(1.0 + eta*eta/(pk1*pk1));
  684.     pk = 2.0 + pk;
  685.   }
  686.  
  687.   D  = 1.0/((pk + pk1)*(x_inv + ek/pk1));
  688.   df = -fcl*(1.0 + ek*ek)*D;
  689.   
  690.   if(fcl != 1.0) fcl = -1.0;
  691.   if(D    < 0.0) fcl = -fcl;
  692.   
  693.   F = F + df;
  694.  
  695.   p = 1.0;
  696.   do {
  697.     pk = pk1;
  698.     pk1 = pk + 1.0;
  699.     ek  = eta / pk;
  700.     tk  = (pk + pk1)*(x_inv + ek/pk1);
  701.     D   =  tk - D*(1.0+ek*ek);
  702.     if(fabs(D) < sqrt(CF1_acc)) {
  703.       p += 1.0;
  704.       if(p > 2.0) {
  705.         printf("HELP............\n");
  706.       }
  707.     }
  708.     D = 1.0/D;
  709.     if(D < 0.0) {
  710.       /* sign of result depends on sign of denominator */
  711.       fcl = -fcl;
  712.     }
  713.     df = df*(D*tk - 1.0);
  714.     F  = F + df;
  715.     if( pk > px ) {
  716.       GSL_ERROR ("error", GSL_ERUNAWAY);
  717.     }
  718.   }
  719.   while(fabs(df) > fabs(F)*CF1_acc);
  720.   
  721.   *fcl_sign = fcl;
  722.   *result = F;
  723.   return GSL_SUCCESS;
  724. }
  725. #endif /* 0 */
  726.  
  727.  
  728. /* Evaluate the second continued fraction to 
  729.  * obtain the ratio
  730.  *    (G' + i F')/(G + i F) := P + i Q
  731.  * at the specified lambda value.
  732.  */
  733. static
  734. int
  735. coulomb_CF2(const double lambda, const double eta, const double x,
  736.             double * result_P, double * result_Q, int * count
  737.             )
  738. {
  739.   int status = GSL_SUCCESS;
  740.  
  741.   const double CF2_acc   = 4.0*GSL_DBL_EPSILON;
  742.   const double CF2_abort = 2.0e+05;
  743.  
  744.   const double wi    = 2.0*eta;
  745.   const double x_inv = 1.0/x;
  746.   const double e2mm1 = eta*eta + lambda*(lambda + 1.0);
  747.   
  748.   double ar = -e2mm1;
  749.   double ai =  eta;
  750.  
  751.   double br =  2.0*(x - eta);
  752.   double bi =  2.0;
  753.  
  754.   double dr =  br/(br*br + bi*bi);
  755.   double di = -bi/(br*br + bi*bi);
  756.  
  757.   double dp = -x_inv*(ar*di + ai*dr);
  758.   double dq =  x_inv*(ar*dr - ai*di);
  759.  
  760.   double A, B, C, D;
  761.  
  762.   double pk =  0.0;
  763.   double P  =  0.0;
  764.   double Q  =  1.0 - eta*x_inv;
  765.  
  766.   *count = 0;
  767.  
  768.   do {
  769.     P += dp;
  770.     Q += dq;
  771.     pk += 2.0;
  772.     ar += pk;
  773.     ai += wi;
  774.     bi += 2.0;
  775.     D  = ar*dr - ai*di + br;
  776.     di = ai*dr + ar*di + bi;
  777.     C  = 1.0/(D*D + di*di);
  778.     dr =  C*D;
  779.     di = -C*di;
  780.     A  = br*dr - bi*di - 1.;
  781.     B  = bi*dr + br*di;
  782.     C  = dp*A  - dq*B;
  783.     dq = dp*B  + dq*A;
  784.     dp = C;
  785.     if(pk > CF2_abort) {
  786.       status = GSL_ERUNAWAY;
  787.       break;
  788.     }
  789.     ++(*count);
  790.   }
  791.   while(fabs(dp)+fabs(dq) > (fabs(P)+fabs(Q))*CF2_acc);
  792.  
  793.   if(Q < CF2_abort*GSL_DBL_EPSILON*fabs(P)) {
  794.     status = GSL_ELOSS;
  795.   }
  796.  
  797.   *result_P = P;
  798.   *result_Q = Q;
  799.   return status;
  800. }
  801.  
  802.  
  803. /* WKB evaluation of F, G. Assumes  0 < x < turning point.
  804.  * Overflows are trapped, GSL_EOVRFLW is signalled,
  805.  * and an exponent is returned such that:
  806.  *
  807.  *   result_F = fjwkb * exp(-exponent)
  808.  *   result_G = gjwkb * exp( exponent)
  809.  *
  810.  * See [Biedenharn et al. Phys. Rev. 97, 542-554 (1955), Section IV]
  811.  *
  812.  * Unfortunately, this is not very accurate in general. The
  813.  * test cases typically have 3-4 digits of precision. One could
  814.  * argue that this is ok for general use because, for instance,
  815.  * F is exponentially small in this region and so the absolute
  816.  * accuracy is still roughly acceptable. But it would be better
  817.  * to have a systematic method for improving the precision. See
  818.  * the Abad+Sesma method discussion below.
  819.  */
  820. static
  821. int
  822. coulomb_jwkb(const double lam, const double eta, const double x,
  823.              gsl_sf_result * fjwkb, gsl_sf_result * gjwkb,
  824.          double * exponent)
  825. {
  826.   const double llp1      = lam*(lam+1.0) + 6.0/35.0;
  827.   const double llp1_eff  = GSL_MAX(llp1, 0.0);
  828.   const double rho_ghalf = sqrt(x*(2.0*eta - x) + llp1_eff);
  829.   const double sinh_arg  = sqrt(llp1_eff/(eta*eta+llp1_eff)) * rho_ghalf / x;
  830.   const double sinh_inv  = log(sinh_arg + sqrt(1.0 + sinh_arg*sinh_arg));
  831.  
  832.   const double phi = fabs(rho_ghalf - eta*atan2(rho_ghalf,x-eta) - sqrt(llp1_eff) * sinh_inv);
  833.  
  834.   const double zeta_half = pow(3.0*phi/2.0, 1.0/3.0);
  835.   const double prefactor = sqrt(M_PI*phi*x/(6.0 * rho_ghalf));
  836.   
  837.   double F = prefactor * 3.0/zeta_half;
  838.   double G = prefactor * 3.0/zeta_half; /* Note the sqrt(3) from Bi normalization */
  839.   double F_exp;
  840.   double G_exp;
  841.   
  842.   const double airy_scale_exp = phi;
  843.   gsl_sf_result ai;
  844.   gsl_sf_result bi;
  845.   gsl_sf_airy_Ai_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &ai);
  846.   gsl_sf_airy_Bi_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &bi);
  847.   F *= ai.val;
  848.   G *= bi.val;
  849.   F_exp = log(F) - airy_scale_exp;
  850.   G_exp = log(G) + airy_scale_exp;
  851.  
  852.   if(G_exp >= GSL_LOG_DBL_MAX) {
  853.     fjwkb->val = F;
  854.     gjwkb->val = G;
  855.     fjwkb->err = 1.0e-3 * fabs(F); /* FIXME: real error here ... could be smaller */
  856.     gjwkb->err = 1.0e-3 * fabs(G);
  857.     *exponent = airy_scale_exp;
  858.     GSL_ERROR ("error", GSL_EOVRFLW);
  859.   }
  860.   else {
  861.     fjwkb->val = exp(F_exp);
  862.     gjwkb->val = exp(G_exp);
  863.     fjwkb->err = 1.0e-3 * fabs(fjwkb->val);
  864.     gjwkb->err = 1.0e-3 * fabs(gjwkb->val);
  865.     *exponent = 0.0;
  866.     return GSL_SUCCESS;
  867.   }
  868. }
  869.  
  870.  
  871. /* Asymptotic evaluation of F and G below the minimal turning point.
  872.  *
  873.  * This is meant to be a drop-in replacement for coulomb_jwkb().
  874.  * It uses the expressions in [Abad+Sesma]. This requires some
  875.  * work because I am not sure where it is valid. They mumble
  876.  * something about |x| < |lam|^(-1/2) or 8|eta x| > lam when |x| < 1.
  877.  * This seems true, but I thought the result was based on a uniform
  878.  * expansion and could be controlled by simply using more terms.
  879.  */
  880. #if 0
  881. static
  882. int
  883. coulomb_AS_xlt2eta(const double lam, const double eta, const double x,
  884.                    gsl_sf_result * f_AS, gsl_sf_result * g_AS,
  885.                double * exponent)
  886. {
  887.   /* no time to do this now... */
  888. }
  889. #endif /* 0 */
  890.  
  891.  
  892.  
  893. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  894.  
  895. int
  896. gsl_sf_coulomb_wave_FG_e(const double eta, const double x,
  897.                             const double lam_F,
  898.                 const int  k_lam_G,      /* lam_G = lam_F - k_lam_G */
  899.                             gsl_sf_result * F, gsl_sf_result * Fp,
  900.                 gsl_sf_result * G, gsl_sf_result * Gp,
  901.                 double * exp_F, double * exp_G)
  902. {
  903.   const double lam_G = lam_F - k_lam_G;
  904.  
  905.   if(x < 0.0 || lam_F <= -0.5 || lam_G <= -0.5) {
  906.     GSL_SF_RESULT_SET(F,  0.0, 0.0);
  907.     GSL_SF_RESULT_SET(Fp, 0.0, 0.0);
  908.     GSL_SF_RESULT_SET(G,  0.0, 0.0);
  909.     GSL_SF_RESULT_SET(Gp, 0.0, 0.0);
  910.     *exp_F = 0.0;
  911.     *exp_G = 0.0;
  912.     GSL_ERROR ("domain error", GSL_EDOM);
  913.   }
  914.   else if(x == 0.0) {
  915.     gsl_sf_result C0;
  916.     CLeta(0.0, eta, &C0);
  917.     GSL_SF_RESULT_SET(F,  0.0, 0.0);
  918.     GSL_SF_RESULT_SET(Fp, 0.0, 0.0);
  919.     GSL_SF_RESULT_SET(G,  0.0, 0.0); /* FIXME: should be Inf */
  920.     GSL_SF_RESULT_SET(Gp, 0.0, 0.0); /* FIXME: should be Inf */
  921.     *exp_F = 0.0;
  922.     *exp_G = 0.0;
  923.     if(lam_F == 0.0){
  924.       GSL_SF_RESULT_SET(Fp, C0.val, C0.err);
  925.     }
  926.     if(lam_G == 0.0) {
  927.       GSL_SF_RESULT_SET(Gp, 1.0/C0.val, fabs(C0.err/C0.val)/fabs(C0.val));
  928.     }
  929.     GSL_ERROR ("domain error", GSL_EDOM);
  930.     /* After all, since we are asking for G, this is a domain error... */
  931.   }
  932.   else if(x < 1.2 && 2.0*M_PI*eta < 0.9*(-GSL_LOG_DBL_MIN) && fabs(eta*x) < 10.0) {
  933.     /* Reduce to a small lambda value and use the series
  934.      * representations for F and G. We cannot allow eta to
  935.      * be large and positive because the connection formula
  936.      * for G_lam is badly behaved due to an underflow in sin(phi_lam) 
  937.      * [see coulomb_FG_series() and coulomb_connection() above].
  938.      * Note that large negative eta is ok however.
  939.      */
  940.     const double SMALL = GSL_SQRT_DBL_EPSILON;
  941.     const int N    = (int)(lam_F + 0.5);
  942.     const int span = GSL_MAX(k_lam_G, N);
  943.     const double lam_min = lam_F - N;    /* -1/2 <= lam_min < 1/2 */
  944.     double F_lam_F, Fp_lam_F;
  945.     double G_lam_G, Gp_lam_G;
  946.     double F_lam_F_err, Fp_lam_F_err;
  947.     double Fp_over_F_lam_F;
  948.     double F_sign_lam_F;
  949.     double F_lam_min_unnorm, Fp_lam_min_unnorm;
  950.     double Fp_over_F_lam_min;
  951.     gsl_sf_result F_lam_min;
  952.     gsl_sf_result G_lam_min, Gp_lam_min;
  953.     double F_scale;
  954.     double Gerr_frac;
  955.     double F_scale_frac_err;
  956.     double F_unnorm_frac_err;
  957.  
  958.     /* Determine F'/F at lam_F. */
  959.     int CF1_count;
  960.     int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
  961.  
  962.     int stat_ser;
  963.     int stat_Fr;
  964.     int stat_Gr;
  965.  
  966.     /* Recurse down with unnormalized F,F' values. */
  967.     F_lam_F  = SMALL;
  968.     Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
  969.     if(span != 0) {
  970.       stat_Fr = coulomb_F_recur(lam_min, span, eta, x,
  971.                                 F_lam_F, Fp_lam_F,
  972.                         &F_lam_min_unnorm, &Fp_lam_min_unnorm
  973.                         );
  974.     }
  975.     else {
  976.       F_lam_min_unnorm  =  F_lam_F;
  977.       Fp_lam_min_unnorm = Fp_lam_F;
  978.       stat_Fr = GSL_SUCCESS;
  979.     }
  980.  
  981.     /* Determine F and G at lam_min. */
  982.     if(lam_min == -0.5) {
  983.       stat_ser = coulomb_FGmhalf_series(eta, x, &F_lam_min, &G_lam_min);
  984.     }
  985.     else if(lam_min == 0.0) {
  986.       stat_ser = coulomb_FG0_series(eta, x, &F_lam_min, &G_lam_min);
  987.     }
  988.     else if(lam_min == 0.5) {
  989.       /* This cannot happen. */
  990.       F->val  = F_lam_F;
  991.       F->err  = 2.0 * GSL_DBL_EPSILON * fabs(F->val);
  992.       Fp->val = Fp_lam_F;
  993.       Fp->err = 2.0 * GSL_DBL_EPSILON * fabs(Fp->val);
  994.       G->val  = G_lam_G;
  995.       G->err  = 2.0 * GSL_DBL_EPSILON * fabs(G->val);
  996.       Gp->val = Gp_lam_G;
  997.       Gp->err = 2.0 * GSL_DBL_EPSILON * fabs(Gp->val);
  998.       *exp_F = 0.0;
  999.       *exp_G = 0.0;
  1000.       GSL_ERROR ("error", GSL_ESANITY);
  1001.     }
  1002.     else {
  1003.       stat_ser = coulomb_FG_series(lam_min, eta, x, &F_lam_min, &G_lam_min);
  1004.     }
  1005.  
  1006.     /* Determine remaining quantities. */
  1007.     Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;
  1008.     Gp_lam_min.val  = Fp_over_F_lam_min*G_lam_min.val - 1.0/F_lam_min.val;
  1009.     Gp_lam_min.err  = fabs(Fp_over_F_lam_min)*G_lam_min.err;
  1010.     Gp_lam_min.err += fabs(1.0/F_lam_min.val) * fabs(F_lam_min.err/F_lam_min.val);
  1011.     F_scale     = F_lam_min.val / F_lam_min_unnorm;
  1012.  
  1013.     /* Apply scale to the original F,F' values. */
  1014.     F_scale_frac_err  = fabs(F_lam_min.err/F_lam_min.val);
  1015.     F_unnorm_frac_err = 2.0*GSL_DBL_EPSILON*(CF1_count+span+1);
  1016.     F_lam_F     *= F_scale;
  1017.     F_lam_F_err  = fabs(F_lam_F) * (F_unnorm_frac_err + F_scale_frac_err);
  1018.     Fp_lam_F    *= F_scale;
  1019.     Fp_lam_F_err = fabs(Fp_lam_F) * (F_unnorm_frac_err + F_scale_frac_err);
  1020.  
  1021.     /* Recurse up to get the required G,G' values. */
  1022.     stat_Gr = coulomb_G_recur(lam_min, GSL_MAX(N-k_lam_G,0), eta, x,
  1023.                               G_lam_min.val, Gp_lam_min.val,
  1024.                       &G_lam_G, &Gp_lam_G
  1025.                       );
  1026.  
  1027.     F->val  = F_lam_F;
  1028.     F->err  = F_lam_F_err;
  1029.     F->err += 2.0 * GSL_DBL_EPSILON * fabs(F_lam_F);
  1030.  
  1031.     Fp->val  = Fp_lam_F;
  1032.     Fp->err  = Fp_lam_F_err;
  1033.     Fp->err += 2.0 * GSL_DBL_EPSILON * fabs(Fp_lam_F);
  1034.  
  1035.     Gerr_frac = fabs(G_lam_min.err/G_lam_min.val) + fabs(Gp_lam_min.err/Gp_lam_min.val);
  1036.  
  1037.     G->val  = G_lam_G;
  1038.     G->err  = Gerr_frac * fabs(G_lam_G);
  1039.     G->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(G->val);
  1040.  
  1041.     Gp->val  = Gp_lam_G;
  1042.     Gp->err  = Gerr_frac * fabs(Gp->val);
  1043.     Gp->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(Gp->val);
  1044.  
  1045.     *exp_F = 0.0;
  1046.     *exp_G = 0.0;
  1047.  
  1048.     return GSL_ERROR_SELECT_4(stat_ser, stat_CF1, stat_Fr, stat_Gr);
  1049.   }
  1050.   else if(x < 2.0*eta) {
  1051.     /* Use WKB approximation to obtain F and G at the two
  1052.      * lambda values, and use the Wronskian and the
  1053.      * continued fractions for F'/F to obtain F' and G'.
  1054.      */
  1055.     gsl_sf_result F_lam_F, G_lam_F;
  1056.     gsl_sf_result F_lam_G, G_lam_G;
  1057.     double exp_lam_F, exp_lam_G;
  1058.     int stat_lam_F;
  1059.     int stat_lam_G;
  1060.     int stat_CF1_lam_F;
  1061.     int stat_CF1_lam_G;
  1062.     int CF1_count;
  1063.     double Fp_over_F_lam_F;
  1064.     double Fp_over_F_lam_G;
  1065.     double F_sign_lam_F;
  1066.     double F_sign_lam_G;
  1067.  
  1068.     stat_lam_F = coulomb_jwkb(lam_F, eta, x, &F_lam_F, &G_lam_F, &exp_lam_F);
  1069.     if(k_lam_G == 0) {
  1070.       stat_lam_G = stat_lam_F;
  1071.       F_lam_G = F_lam_F;
  1072.       G_lam_G = G_lam_F;
  1073.       exp_lam_G = exp_lam_F;
  1074.     }
  1075.     else {
  1076.       stat_lam_G = coulomb_jwkb(lam_G, eta, x, &F_lam_G, &G_lam_G, &exp_lam_G);
  1077.     }
  1078.  
  1079.     stat_CF1_lam_F = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
  1080.     if(k_lam_G == 0) {
  1081.       stat_CF1_lam_G  = stat_CF1_lam_F;
  1082.       F_sign_lam_G    = F_sign_lam_F;
  1083.       Fp_over_F_lam_G = Fp_over_F_lam_F;
  1084.     }
  1085.     else {
  1086.       stat_CF1_lam_G = coulomb_CF1(lam_G, eta, x, &F_sign_lam_G, &Fp_over_F_lam_G, &CF1_count);
  1087.     }
  1088.  
  1089.     F->val = F_lam_F.val;
  1090.     F->err = F_lam_F.err;
  1091.  
  1092.     G->val = G_lam_G.val;
  1093.     G->err = G_lam_G.err;
  1094.  
  1095.     Fp->val  = Fp_over_F_lam_F * F_lam_F.val;
  1096.     Fp->err  = fabs(Fp_over_F_lam_F) * F_lam_F.err;
  1097.     Fp->err += 2.0*GSL_DBL_EPSILON*fabs(Fp->val);
  1098.  
  1099.     Gp->val  = Fp_over_F_lam_G * G_lam_G.val - 1.0/F_lam_G.val;
  1100.     Gp->err  = fabs(Fp_over_F_lam_G) * G_lam_G.err;
  1101.     Gp->err += fabs(1.0/F_lam_G.val) * fabs(F_lam_G.err/F_lam_G.val);
  1102.  
  1103.     *exp_F = exp_lam_F;
  1104.     *exp_G = exp_lam_G;
  1105.  
  1106.     if(stat_lam_F == GSL_EOVRFLW || stat_lam_G == GSL_EOVRFLW) {
  1107.       GSL_ERROR ("overflow", GSL_EOVRFLW);
  1108.     }
  1109.     else {
  1110.       return GSL_ERROR_SELECT_2(stat_lam_F, stat_lam_G);
  1111.     }
  1112.   }
  1113.   else {
  1114.     /* x > 2 eta, so we know that we can find a lambda value such
  1115.      * that x is above the turning point. We do this, evaluate
  1116.      * using Steed's method at that oscillatory point, then
  1117.      * use recursion on F and G to obtain the required values.
  1118.      *
  1119.      * lam_0   = a value of lambda such that x is below the turning point
  1120.      * lam_min = minimum of lam_0 and the requested lam_G, since
  1121.      *           we must go at least as low as lam_G
  1122.      */
  1123.     const double SMALL = GSL_SQRT_DBL_EPSILON;
  1124.     const double C = sqrt(1.0 + 4.0*x*(x-2.0*eta));
  1125.     const int N = ceil(lam_F - C + 0.5);
  1126.     const double lam_0   = lam_F - GSL_MAX(N, 0);
  1127.     const double lam_min = GSL_MIN(lam_0, lam_G);
  1128.     double F_lam_F, Fp_lam_F;
  1129.     double G_lam_G, Gp_lam_G;
  1130.     double F_lam_min_unnorm, Fp_lam_min_unnorm;
  1131.     double F_lam_min, Fp_lam_min;
  1132.     double G_lam_min, Gp_lam_min;
  1133.     double Fp_over_F_lam_F;
  1134.     double Fp_over_F_lam_min;
  1135.     double F_sign_lam_F;
  1136.     double P_lam_min, Q_lam_min;
  1137.     double alpha;
  1138.     double gamma;
  1139.     double F_scale;
  1140.  
  1141.     int CF1_count;
  1142.     int CF2_count;
  1143.     int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
  1144.     int stat_CF2;
  1145.     int stat_Fr;
  1146.     int stat_Gr;
  1147.  
  1148.     int F_recur_count;
  1149.     int G_recur_count;
  1150.  
  1151.     double err_amplify;
  1152.  
  1153.     F_lam_F  = SMALL;
  1154.     Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
  1155.  
  1156.     /* Backward recurrence to get F,Fp at lam_min */
  1157.     F_recur_count = GSL_MAX(k_lam_G, N);
  1158.     stat_Fr = coulomb_F_recur(lam_min, F_recur_count, eta, x,
  1159.                               F_lam_F, Fp_lam_F,
  1160.                       &F_lam_min_unnorm, &Fp_lam_min_unnorm
  1161.                       );
  1162.     Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;
  1163.  
  1164.     /* Steed evaluation to complete evaluation of F,Fp,G,Gp at lam_min */
  1165.     stat_CF2 = coulomb_CF2(lam_min, eta, x, &P_lam_min, &Q_lam_min, &CF2_count);
  1166.     alpha = Fp_over_F_lam_min - P_lam_min;
  1167.     gamma = alpha/Q_lam_min;
  1168.     F_lam_min  = F_sign_lam_F / sqrt(alpha*alpha/Q_lam_min + Q_lam_min);
  1169.     Fp_lam_min = Fp_over_F_lam_min * F_lam_min;
  1170.     G_lam_min  = gamma * F_lam_min;
  1171.     Gp_lam_min = (P_lam_min * gamma - Q_lam_min) * F_lam_min;
  1172.  
  1173.     /* Apply scale to values of F,Fp at lam_F (the top). */
  1174.     F_scale = F_lam_min / F_lam_min_unnorm;    
  1175.     F_lam_F  *= F_scale;
  1176.     Fp_lam_F *= F_scale;
  1177.  
  1178.     /* Forward recurrence to get G,Gp at lam_G (the top). */
  1179.     G_recur_count = GSL_MAX(N-k_lam_G,0);
  1180.     stat_Gr = coulomb_G_recur(lam_min, G_recur_count, eta, x,
  1181.                               G_lam_min, Gp_lam_min,
  1182.                       &G_lam_G, &Gp_lam_G
  1183.                       );
  1184.  
  1185.     err_amplify = CF1_count + CF2_count + F_recur_count + G_recur_count + 1;
  1186.  
  1187.     F->val  = F_lam_F;
  1188.     F->err  = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(F->val);
  1189.  
  1190.     Fp->val = Fp_lam_F;
  1191.     Fp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Fp->val);
  1192.  
  1193.     G->val  = G_lam_G;
  1194.     G->err  = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(G->val);
  1195.  
  1196.     Gp->val = Gp_lam_G;
  1197.     Gp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Gp->val);
  1198.  
  1199.     *exp_F = 0.0;
  1200.     *exp_G = 0.0;
  1201.  
  1202.     return GSL_ERROR_SELECT_4(stat_CF1, stat_CF2, stat_Fr, stat_Gr);
  1203.   }
  1204. }
  1205.  
  1206.  
  1207. int
  1208. gsl_sf_coulomb_wave_F_array(double lam_min, int kmax,
  1209.                                  double eta, double x, 
  1210.                                  double * fc_array,
  1211.                      double * F_exp)
  1212. {
  1213.   if(x == 0.0) {
  1214.     int k;
  1215.     *F_exp = 0.0;
  1216.     for(k=0; k<=kmax; k++) {
  1217.       fc_array[k] = 0.0;
  1218.     }
  1219.     if(lam_min == 0.0){
  1220.       gsl_sf_result f_0;
  1221.       CLeta(0.0, eta, &f_0);
  1222.       fc_array[0] = f_0.val;
  1223.     }
  1224.     return GSL_SUCCESS;
  1225.   }
  1226.   else {
  1227.     const double x_inv = 1.0/x;
  1228.     const double lam_max = lam_min + kmax;
  1229.     gsl_sf_result F, Fp;
  1230.     gsl_sf_result G, Gp;
  1231.     double G_exp;
  1232.  
  1233.     int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, 0,
  1234.                                               &F, &Fp, &G, &Gp, F_exp, &G_exp);
  1235.  
  1236.     double fcl  = F.val;
  1237.     double fpl = Fp.val;
  1238.     double lam = lam_max;
  1239.     int k;
  1240.  
  1241.     fc_array[kmax] = F.val;
  1242.  
  1243.     for(k=kmax-1; k>=0; k--) {
  1244.       double el = eta/lam;
  1245.       double rl = sqrt(1.0 + el*el);
  1246.       double sl = el  + lam*x_inv;
  1247.       double fc_lm1 = (fcl*sl + fpl)/rl;
  1248.       fc_array[k]   = fc_lm1;
  1249.       fpl           =  fc_lm1*sl - fcl*rl;
  1250.       fcl           =  fc_lm1;
  1251.       lam -= 1.0;
  1252.     }
  1253.  
  1254.     return stat_FG;
  1255.   }
  1256. }
  1257.  
  1258.  
  1259. int
  1260. gsl_sf_coulomb_wave_FG_array(double lam_min, int kmax,
  1261.                                   double eta, double x,
  1262.                                   double * fc_array, double * gc_array,
  1263.                       double * F_exp, double * G_exp)
  1264. {
  1265.   const double x_inv = 1.0/x;
  1266.   const double lam_max = lam_min + kmax;
  1267.   gsl_sf_result F, Fp;
  1268.   gsl_sf_result G, Gp;
  1269.  
  1270.   int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
  1271.                                             &F, &Fp, &G, &Gp, F_exp, G_exp);
  1272.  
  1273.   double fcl  = F.val;
  1274.   double fpl = Fp.val;
  1275.   double lam = lam_max;
  1276.   int k;
  1277.  
  1278.   double gcl, gpl;
  1279.  
  1280.   fc_array[kmax] = F.val;
  1281.  
  1282.   for(k=kmax-1; k>=0; k--) {
  1283.     double el = eta/lam;
  1284.     double rl = sqrt(1.0 + el*el);
  1285.     double sl = el  + lam*x_inv;
  1286.     double fc_lm1;
  1287.     fc_lm1 = (fcl*sl + fpl)/rl;
  1288.     fc_array[k] = fc_lm1;
  1289.     fpl         =  fc_lm1*sl - fcl*rl;
  1290.     fcl         =  fc_lm1;
  1291.     lam -= 1.0;
  1292.   }
  1293.  
  1294.   gcl = G.val;
  1295.   gpl = Gp.val;
  1296.   lam = lam_min + 1.0;
  1297.  
  1298.   gc_array[0] = G.val;
  1299.  
  1300.   for(k=1; k<=kmax; k++) {
  1301.     double el = eta/lam;
  1302.     double rl = sqrt(1.0 + el*el);
  1303.     double sl = el + lam*x_inv;
  1304.     double gcl1 = (sl*gcl - gpl)/rl;
  1305.     gc_array[k] = gcl1;
  1306.     gpl         = rl*gcl - sl*gcl1;
  1307.     gcl         = gcl1;
  1308.     lam += 1.0;
  1309.   }
  1310.  
  1311.   return stat_FG;
  1312. }
  1313.  
  1314.  
  1315. int
  1316. gsl_sf_coulomb_wave_FGp_array(double lam_min, int kmax,
  1317.                                    double eta, double x,
  1318.                            double * fc_array, double * fcp_array,
  1319.                            double * gc_array, double * gcp_array,
  1320.                        double * F_exp, double * G_exp)
  1321.  
  1322. {
  1323.   const double x_inv = 1.0/x;
  1324.   const double lam_max = lam_min + kmax;
  1325.   gsl_sf_result F, Fp;
  1326.   gsl_sf_result G, Gp;
  1327.  
  1328.   int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
  1329.                                             &F, &Fp, &G, &Gp, F_exp, G_exp);
  1330.  
  1331.   double fcl  = F.val;
  1332.   double fpl = Fp.val;
  1333.   double lam = lam_max;
  1334.   int k;
  1335.  
  1336.   double gcl, gpl;
  1337.  
  1338.   fc_array[kmax]  = F.val;
  1339.   fcp_array[kmax] = Fp.val;
  1340.  
  1341.   for(k=kmax-1; k>=0; k--) {
  1342.     double el = eta/lam;
  1343.     double rl = sqrt(1.0 + el*el);
  1344.     double sl = el  + lam*x_inv;
  1345.     double fc_lm1;
  1346.     fc_lm1 = (fcl*sl + fpl)/rl;
  1347.     fc_array[k]  = fc_lm1;
  1348.     fpl          = fc_lm1*sl - fcl*rl;
  1349.     fcp_array[k] = fpl;
  1350.     fcl          =  fc_lm1;
  1351.     lam -= 1.0;
  1352.   }
  1353.  
  1354.   gcl = G.val;
  1355.   gpl = Gp.val;
  1356.   lam = lam_min + 1.0;
  1357.  
  1358.   gc_array[0]  = G.val;
  1359.   gcp_array[0] = Gp.val;
  1360.  
  1361.   for(k=1; k<=kmax; k++) {
  1362.     double el = eta/lam;
  1363.     double rl = sqrt(1.0 + el*el);
  1364.     double sl = el + lam*x_inv;
  1365.     double gcl1 = (sl*gcl - gpl)/rl;
  1366.     gc_array[k]  = gcl1;
  1367.     gpl          = rl*gcl - sl*gcl1;
  1368.     gcp_array[k] = gpl;
  1369.     gcl          = gcl1;
  1370.     lam += 1.0;
  1371.   }
  1372.  
  1373.   return stat_FG;
  1374. }
  1375.  
  1376.  
  1377. int
  1378. gsl_sf_coulomb_wave_sphF_array(double lam_min, int kmax,
  1379.                                     double eta, double x,
  1380.                             double * fc_array,
  1381.                         double * F_exp)
  1382. {
  1383.   int k;
  1384.  
  1385.   if(x < 0.0 || lam_min < -0.5) {
  1386.     GSL_ERROR ("domain error", GSL_EDOM);
  1387.   }
  1388.   else if(x < 10.0/GSL_DBL_MAX) {
  1389.     for(k=0; k<=kmax; k++) {
  1390.       fc_array[k] = 0.0;
  1391.     }
  1392.     if(lam_min == 0.0) {
  1393.       fc_array[0] = sqrt(C0sq(eta));
  1394.     }
  1395.     *F_exp = 0.0;
  1396.     if(x == 0.0)
  1397.       return GSL_SUCCESS;
  1398.     else
  1399.       GSL_ERROR ("underflow", GSL_EUNDRFLW);
  1400.   }
  1401.   else {
  1402.     int k;
  1403.     int stat_F = gsl_sf_coulomb_wave_F_array(lam_min, kmax,
  1404.                                                   eta, x, 
  1405.                                                   fc_array,
  1406.                                                   F_exp);
  1407.  
  1408.     for(k=0; k<=kmax; k++) {
  1409.       fc_array[k] = fc_array[k] / x;
  1410.     }
  1411.     return stat_F;
  1412.   }
  1413. }
  1414.  
  1415.  
  1416.